Integrating Human Factors Expertise into Development of Automated Vehicles
Replication package
Amna Pir Muhammad, Alessia Knauss, Eric Knauss, Jonas Bärgman
First revision 2024-04-25. This revision 2024-04-29.
1 Introduction
This document serves as a comprehensive guide to replicate the analysis conducted in our experiments. The replication package enables anyone to download, rerun the analysis, and verify our assumptions.
1.1 Data Loading
To begin, load the dataset available in the Data/ directory of our GitHub repository. The dataset is read into R using the following code:
data <- read_sav("Untitled26.sav")In our dataset: Ranking Variable: Presented by the M201_ prefix, this variable assesses preferences across ten options. Effectiveness Variable: Presented by the M209_ prefix, this variable measures the effectiveness of different options. Easiness Variable: Presented by the M210_ prefix, this variable evaluates the perceived ease of each option.
For more detailed information on these variables and their corresponding options, please refer to our manuscript or questionnaire.
2 Model Design
2.1 Effectiveness
2.1.1 Convert Categorical Variables to Factors
We convert specific columns to ordered factors to ensure correct ordering in categorical analysis.
data$D104 <- factor(data$D104, levels = c("1", "2", "3"), ordered = TRUE)
data$D102 <- factor(data$D102, levels = c("1", "2", "3", "4"), ordered = TRUE)2.1.2 Data Cleaning - Handle Missing Values
We handle missing values by recoding -9 to NA for the variables starting with “M201_,”M209_“, and”M210_“.
modeldata <- data %>%
mutate(across(starts_with("M201_"), ~if_else(. == -9, NA_real_, .))) %>%
mutate(across(starts_with("M209_"), ~if_else(. == -9, NA_real_, .))) %>%
mutate(across(starts_with("M210_"), ~if_else(. == -9, NA_real_, .)))2.1.3 Create New Variables Based on Existing Ones
We create new variables based on existing ones to simplify and categorize the data for analysis.
modeldata <- modeldata %>%
mutate(
D104 = D104,
Expertise = case_when(
D102 == 1 ~ "Human Factors",
D102 == 2 ~ "Engineering",
D102 == 3 ~ "Both",
TRUE ~ "Other"
),
Dev_Struc = case_when(
I004_01 >= 1 & I004_01 <= 35 ~ "Waterfall",
I004_01 >= 36 & I004_01 <= 70 ~ "Hybrid",
I004_01 >= 71 & I004_01 <= 100 ~ "Agile",
TRUE ~ "NA"
),
Experience = case_when(
D104 == 1 ~ "LevelLow",
D104 == 2 ~ "LevelMedium",
D104 == 3 ~ "LevelHigh",
TRUE ~ "Other"
),
AreaOfWork = case_when(
D101_01 == 2 ~ "Supplier",
D101_02 == 2 ~ "OEM",
D101_03 == 2 ~ "Consultancy",
D101_04 == 2 ~ "Research Institute",
D101_05 == 2 ~ "University",
D101_06 == 2 ~ "Government",
D101_07 == 2 ~ "Non-Automotive",
TRUE ~ "Other"
)
) %>%
pivot_longer(cols = starts_with("M209_"), names_to = "placementOption", values_to = "effectivenessRating")2.1.4 Convert Effectiveness varaible to Factor
We convert the “effectivenessRating” variable to an ordered factor for correct ordering in the analysis.
modeldata$effectivenessRating <- factor(modeldata$effectivenessRating, levels = c("1", "2", "3", "4", "5"), ordered = TRUE)2.1.5 Define Formula
We define the formula to predict “Effectiveness” based on various predictors.
model_formula <- bf(effectivenessRating ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1|ID))2.1.6 Prior Specification
We specify the prior distributions for the model parameters.
prior_spec <- c(
prior(normal(0, 5), "b"),
prior(normal(0, 10), "Intercept"),
prior(cauchy(0, 2), "sd")
)2.1.7 Prior Predictive Checks
We check the and plot the priors from the sample only without empirical data. Sample only from the priors. Evidently the medians are quite evenly set along the \(x\)-axis, and the uncertainty is fairly uniformly distributed among the categories \(1,\ldots,6\) (the bars).
fit <- brm(
formula = model_formula,
data = modeldata,
family = cumulative(probit),
prior = prior_spec,
chains = 4,
sample_prior = "only",
cores = 4,
seed = 12345
)## Warning: Rows containing NAs were excluded from the model.## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored## Start sampling## Warning: 1 of 4000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.pp_check(fit, type = "bars", nsamples = 250)## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.
### Sample with data
We fit the model to the data using the specified formula, priors, and settings.
fit <- brm(
formula = model_formula,
data = modeldata,
family = cumulative(probit),
prior = prior_spec,
chains = 4,
cores = 4,
seed = 12345
)## Warning: Rows containing NAs were excluded from the model.## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored## Start sampling2.1.8 Summary and Diagnostics
We print the summary of the fitted model and perform diagnostic checks to ensure the model has converged and is reliable. Our caterpillar plots look good for all estimated parameters (i.e., they look like fat caterpillars when the chains have mixed well).
mcmc_trace(fit, regex_pars = "^b_") + legend_none()
#Diagnostics such as divergences, tree depth, energy, $\widehat{R}$, and $\mathrm{ESS}$ all look good.
rstan::check_hmc_diagnostics(eval(fit)$fit)##
## Divergences:## 0 of 4000 iterations ended with a divergence.##
## Tree depth:## 0 of 4000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.if(max(rhat(eval(fit)), na.rm=T) >= 1.01) {
print("Warning: Rhat >=1.01")
} else {
print("All Rhat <1.01")
}## [1] "All Rhat <1.01"if(min(neff_ratio(eval(fit)), na.rm=T) <= 0.2) {
print("Warning: ESS <=0.2")
} else {
print("All ESS >0.2")
}## [1] "All ESS >0.2"2.1.9 Posterior Predictive Checks
We perform posterior predictive checks to assess how well the model fits the data.
pp_check(fit, type = "bars", nsamples = 250)## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.
2.1.10 Ploting
plot(fit)
## Extract and Plot Posterior Distributions
We extract posterior samples and plot the posterior distributions for the model coefficients.
posterior_samples <- posterior_samples(fit)## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.coefficients <- grep("b_", colnames(posterior_samples), value = TRUE)
plots <- lapply(coefficients, function(coef_name) {
ggplot(posterior_samples, aes_string(x = coef_name)) +
geom_density(fill = "blue", alpha = 0.5) +
labs(title = paste("Posterior Distribution for", gsub("b_", "", coef_name)),
x = "Effect Size", y = "Density")
})## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.#for (p in plots) {
# print(p)
#}2.1.11 Ploting the posterior probability
Next, let’s plot the posterior probability densities for our population-level estimates on the \(\logit\) scale.
parameter_names <- colnames(posterior_samples)
parameters <- grep("^b_", parameter_names, value = TRUE)
parameters <- parameters[!grepl("Intercept", parameters)]
print(parameters)## [1] "b_placementOptionM209_02" "b_placementOptionM209_03"
## [3] "b_placementOptionM209_04" "b_placementOptionM209_05"
## [5] "b_placementOptionM209_06" "b_placementOptionM209_07"
## [7] "b_placementOptionM209_08" "b_placementOptionM209_09"
## [9] "b_placementOptionM209_10" "b_ExperienceLevelLow"
## [11] "b_ExperienceLevelMedium" "b_ExpertiseEngineering"
## [13] "b_ExpertiseHumanFactors" "b_ExpertiseOther"
## [15] "b_AreaOfWorkNonMAutomotive" "b_AreaOfWorkOEM"
## [17] "b_AreaOfWorkResearchInstitute" "b_AreaOfWorkSupplier"
## [19] "b_Dev_StrucHybrid" "b_Dev_StrucWaterfall"parameters <- rev(parameters)
p <- mcmc_areas(fit, pars = parameters, prob_outer = 0.95, prob = 0.5)
#+ vline_0(size = 0.3)
y_labels <- c("Feature Definition", "Feature Realization", "System Design",
"Dedicated human factors team", "User-oriented teams", "Non-functional requirements team",
"System/feature evaluation team", "Safety team", "Person/team responsible for the overall system",
"ExperienceLevelLow", "ExperienceLevelMedium", "ExpertiseEngineering",
"ExpertiseHumanFactors", "ExpertiseOther", "AreaOfWorkNonMAutomotive",
"AreaOfWorkOEM", "AreaOfWorkResearchInstitute", "AreaOfWorkSupplier")
y_labels <- rev(y_labels)
p <- p + scale_y_discrete(labels = y_labels) ## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.print(p)
2.1.12 Conditional effects
Below we plot conditional effects for some of our estimates.
ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "placementOption", categorical = TRUE)
p_crit <- plot(p)[[1]] + # Removed plot = FALSE
scale_color_viridis(discrete = TRUE) +
scale_fill_viridis(discrete = TRUE) +
xlab("placementOption") +
ylab("") +
scale_x_continuous(breaks=seq(0,1))
ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Experience", categorical = TRUE)
p_crit <- plot(p)[[1]] + # Removed plot = FALSE
scale_color_viridis(discrete = TRUE) +
scale_fill_viridis(discrete = TRUE) +
xlab("Experience") +
ylab("") +
scale_x_continuous(breaks=seq(0,1))
ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Expertise", categorical = TRUE)
p_crit <- plot(p)[[1]] + # Removed plot = FALSE
scale_color_viridis(discrete = TRUE) +
scale_fill_viridis(discrete = TRUE) +
xlab("Criticality") +
ylab("") +
scale_x_continuous(breaks=seq(0,1))
2.2 Easiness
2.2.1 Convert Categorical Variables to Factors
We convert specific columns to ordered factors to ensure correct ordering in categorical analysis.
data$D104 <- factor(data$D104, levels = c("1", "2", "3"), ordered = TRUE)
data$D102 <- factor(data$D102, levels = c("1", "2", "3", "4"), ordered = TRUE)2.2.2 Data Cleaning - Handle Missing Values
We handle missing values by recoding -9 to NA for the variables starting with “M201_,”M209_“, and”M210_“.
modeldata <- data %>%
mutate(across(starts_with("M201_"), ~if_else(. == -9, NA_real_, .))) %>%
mutate(across(starts_with("M209_"), ~if_else(. == -9, NA_real_, .))) %>%
mutate(across(starts_with("M210_"), ~if_else(. == -9, NA_real_, .)))2.2.3 Create New Variables Based on Existing Ones
We create new variables based on existing ones to simplify and categorize the data for analysis.
modeldata <- modeldata %>%
mutate(
D104 = D104,
Expertise = case_when(
D102 == 1 ~ "Human Factors",
D102 == 2 ~ "Engineering",
D102 == 3 ~ "Both",
TRUE ~ "Other"
),
Dev_Struc = case_when(
I004_01 >= 1 & I004_01 <= 35 ~ "Waterfall",
I004_01 >= 36 & I004_01 <= 70 ~ "Hybrid",
I004_01 >= 71 & I004_01 <= 100 ~ "Agile",
TRUE ~ "NA"
),
Experience = case_when(
D104 == 1 ~ "LevelLow",
D104 == 2 ~ "LevelMedium",
D104 == 3 ~ "LevelHigh",
TRUE ~ "Other"
),
AreaOfWork = case_when(
D101_01 == 2 ~ "Supplier",
D101_02 == 2 ~ "OEM",
D101_03 == 2 ~ "Consultancy",
D101_04 == 2 ~ "Research Institute",
D101_05 == 2 ~ "University",
D101_06 == 2 ~ "Government",
D101_07 == 2 ~ "Non-Automotive",
TRUE ~ "Other"
)
) %>%
pivot_longer(cols = starts_with("M210_"), names_to = "placementOption", values_to = "easinessRating")2.2.4 Convert Easiness varaible to Factor
We convert the “easinessRating” variable to an ordered factor for correct ordering in the analysis.
modeldata$easinessRating <- factor(modeldata$easinessRating, levels = c("1", "2", "3", "4", "5"), ordered = TRUE)2.2.5 Define Formula
We define the formula to predict “Easiness” based on various predictors.
model_formula <- bf(easinessRating ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1|ID))2.2.6 Prior Specification
We specify the prior distributions for the model parameters.
prior_spec <- c(
prior(normal(0, 5), "b"),
prior(normal(0, 10), "Intercept"),
prior(cauchy(0, 2), "sd")
)2.2.7 Prior Predictive Checks
We check the and plot the priors from the sample only without empirical data. Sample only from the priors. Evidently the medians are quite evenly set along the \(x\)-axis, and the uncertainty is fairly uniformly distributed among the categories \(1,\ldots,6\) (the bars).
fit <- brm(
formula = model_formula,
data = modeldata,
family = cumulative(probit),
prior = prior_spec,
chains = 4,
sample_prior = "only",
cores = 4,
seed = 12345
)## Warning: Rows containing NAs were excluded from the model.## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored## Start sampling## Warning: 1 of 4000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.pp_check(fit, type = "bars", nsamples = 250)## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.
### Sample with data
We fit the model to the data using the specified formula, priors, and settings.
fit <- brm(
formula = model_formula,
data = modeldata,
family = cumulative(probit),
prior = prior_spec,
chains = 4,
cores = 4,
seed = 12345
)## Warning: Rows containing NAs were excluded from the model.## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored## Start sampling2.2.8 Summary and Diagnostics
We print the summary of the fitted model and perform diagnostic checks to ensure the model has converged and is reliable. Our caterpillar plots look good for all estimated parameters (i.e., they look like fat caterpillars when the chains have mixed well).
mcmc_trace(fit, regex_pars = "^b_") + legend_none()
#Diagnostics such as divergences, tree depth, energy, $\widehat{R}$, and $\mathrm{ESS}$ all look good.
rstan::check_hmc_diagnostics(eval(fit)$fit)##
## Divergences:## 0 of 4000 iterations ended with a divergence.##
## Tree depth:## 0 of 4000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.if(max(rhat(eval(fit)), na.rm=T) >= 1.01) {
print("Warning: Rhat >=1.01")
} else {
print("All Rhat <1.01")
}## [1] "All Rhat <1.01"if(min(neff_ratio(eval(fit)), na.rm=T) <= 0.2) {
print("Warning: ESS <=0.2")
} else {
print("All ESS >0.2")
}## [1] "Warning: ESS <=0.2"2.2.9 Posterior Predictive Checks
We perform posterior predictive checks to assess how well the model fits the data.
pp_check(fit, type = "bars", nsamples = 250)## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.
2.2.10 Ploting
plot(fit)
## Extract and Plot Posterior Distributions
We extract posterior samples and plot the posterior distributions for the model coefficients.
posterior_samples <- posterior_samples(fit)## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.coefficients <- grep("b_", colnames(posterior_samples), value = TRUE)
plots <- lapply(coefficients, function(coef_name) {
ggplot(posterior_samples, aes_string(x = coef_name)) +
geom_density(fill = "blue", alpha = 0.5) +
labs(title = paste("Posterior Distribution for", gsub("b_", "", coef_name)),
x = "Effect Size", y = "Density")
})
#for (p in plots) {
# print(p)
#}2.2.11 Ploting the posterior probability
Next, let’s plot the posterior probability densities for our population-level estimates on the \(\logit\) scale.
parameter_names <- colnames(posterior_samples)
parameters <- grep("^b_", parameter_names, value = TRUE)
parameters <- parameters[!grepl("Intercept", parameters)]
print(parameters)## [1] "b_placementOptionM210_02" "b_placementOptionM210_03"
## [3] "b_placementOptionM210_04" "b_placementOptionM210_05"
## [5] "b_placementOptionM210_06" "b_placementOptionM210_07"
## [7] "b_placementOptionM210_08" "b_placementOptionM210_09"
## [9] "b_placementOptionM210_10" "b_ExperienceLevelLow"
## [11] "b_ExperienceLevelMedium" "b_ExpertiseEngineering"
## [13] "b_ExpertiseHumanFactors" "b_ExpertiseOther"
## [15] "b_AreaOfWorkNonMAutomotive" "b_AreaOfWorkOEM"
## [17] "b_AreaOfWorkResearchInstitute" "b_AreaOfWorkSupplier"
## [19] "b_Dev_StrucHybrid" "b_Dev_StrucWaterfall"parameters <- rev(parameters)
p <- mcmc_areas(fit, pars = parameters, prob_outer = 0.95, prob = 0.5)
#+ vline_0(size = 0.3)
y_labels <- c("Feature Definition", "Feature Realization", "System Design",
"Dedicated human factors team", "User-oriented teams", "Non-functional requirements team",
"System/feature evaluation team", "Safety team", "Person/team responsible for the overall system",
"ExperienceLevelLow", "ExperienceLevelMedium", "ExpertiseEngineering",
"ExpertiseHumanFactors", "ExpertiseOther", "AreaOfWorkNonMAutomotive",
"AreaOfWorkOEM", "AreaOfWorkResearchInstitute", "AreaOfWorkSupplier")
y_labels <- rev(y_labels)
p <- p + scale_y_discrete(labels = y_labels) ## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.print(p)
2.2.12 Conditional effects
Below we plot conditional effects for some of our estimates.
ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "placementOption", categorical = TRUE)
p_crit <- plot(p)[[1]] + # Removed plot = FALSE
scale_color_viridis(discrete = TRUE) +
scale_fill_viridis(discrete = TRUE) +
xlab("placementOption") +
ylab("") +
scale_x_continuous(breaks=seq(0,1))
ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Experience", categorical = TRUE)
p_crit <- plot(p)[[1]] + # Removed plot = FALSE
scale_color_viridis(discrete = TRUE) +
scale_fill_viridis(discrete = TRUE) +
xlab("Experience") +
ylab("") +
scale_x_continuous(breaks=seq(0,1))
ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Expertise", categorical = TRUE)
p_crit <- plot(p)[[1]] + # Removed plot = FALSE
scale_color_viridis(discrete = TRUE) +
scale_fill_viridis(discrete = TRUE) +
xlab("Criticality") +
ylab("") +
scale_x_continuous(breaks=seq(0,1))
2.3 Ranking
2.3.1 Convert Categorical Variables to Factors
We convert specific columns to ordered factors to ensure correct ordering in categorical analysis.
data$D104 <- factor(data$D104, levels = c("1", "2", "3"), ordered = TRUE)
data$D102 <- factor(data$D102, levels = c("1", "2", "3", "4"), ordered = TRUE)2.3.2 Data Cleaning - Handle Missing Values
We handle missing values by recoding -9 to NA for the variables starting with “M201_,”M209_“, and”M210_“.
modeldata <- data %>%
mutate(across(starts_with("M201_"), ~if_else(. == -9, NA_real_, .))) %>%
mutate(across(starts_with("M209_"), ~if_else(. == -9, NA_real_, .))) %>%
mutate(across(starts_with("M210_"), ~if_else(. == -9, NA_real_, .)))2.3.3 Create New Variables Based on Existing Ones
We create new variables based on existing ones to simplify and categorize the data for analysis.
modeldata <- modeldata %>%
mutate(
D104 = D104,
Expertise = case_when(
D102 == 1 ~ "Human Factors",
D102 == 2 ~ "Engineering",
D102 == 3 ~ "Both",
TRUE ~ "Other"
),
Dev_Struc = case_when(
I004_01 >= 1 & I004_01 <= 35 ~ "Waterfall",
I004_01 >= 36 & I004_01 <= 70 ~ "Hybrid",
I004_01 >= 71 & I004_01 <= 100 ~ "Agile",
TRUE ~ "NA"
),
Experience = case_when(
D104 == 1 ~ "LevelLow",
D104 == 2 ~ "LevelMedium",
D104 == 3 ~ "LevelHigh",
TRUE ~ "Other"
),
AreaOfWork = case_when(
D101_01 == 2 ~ "Supplier",
D101_02 == 2 ~ "OEM",
D101_03 == 2 ~ "Consultancy",
D101_04 == 2 ~ "Research Institute",
D101_05 == 2 ~ "University",
D101_06 == 2 ~ "Government",
D101_07 == 2 ~ "Non-Automotive",
TRUE ~ "Other"
)
) %>%
pivot_longer(cols = starts_with("M201_"), names_to = "placementOption", values_to = "ranking")2.3.4 Convert Ranking varaible to Factor
We convert the “ranking” variable to an ordered factor for correct ordering in the analysis.
modeldata$ranking <- factor(modeldata$ranking, levels = c("1", "2", "3", "4", "5", "6","7", "8","9", "10"), ordered = TRUE)2.3.5 Define Formula
We define the formula to predict “Ranking” based on various predictors.
model_formula <- bf(ranking ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1|ID))2.3.6 Prior Specification
We specify the prior distributions for the model parameters.
prior_spec <- c(
prior(normal(0, 5), "b"),
prior(normal(0, 10), "Intercept"),
prior(cauchy(0, 2), "sd")
)2.3.7 Prior Predictive Checks
We check the and plot the priors from the sample only without empirical data. Sample only from the priors. Evidently the medians are quite evenly set along the \(x\)-axis, and the uncertainty is fairly uniformly distributed among the categories \(1,\ldots,6\) (the bars).
fit <- brm(
formula = model_formula,
data = modeldata,
family = cumulative(probit),
prior = prior_spec,
chains = 4,
sample_prior = "only",
cores = 4,
seed = 12345
)## Warning: Rows containing NAs were excluded from the model.## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored## Start sampling## Warning: 3 of 4000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.pp_check(fit, type = "bars", nsamples = 250)## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.
### Sample with data
We fit the model to the data using the specified formula, priors, and settings.
fit <- brm(
formula = model_formula,
data = modeldata,
family = cumulative(probit),
prior = prior_spec,
chains = 4,
cores = 4,
seed = 12345
)## Warning: Rows containing NAs were excluded from the model.## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored## Start sampling## Chain 1 Initialization failed.## Chain 4 Initialization failed.2.3.8 Summary and Diagnostics
We print the summary of the fitted model and perform diagnostic checks to ensure the model has converged and is reliable. Our caterpillar plots look good for all estimated parameters (i.e., they look like fat caterpillars when the chains have mixed well).
mcmc_trace(fit, regex_pars = "^b_") + legend_none()
#Diagnostics such as divergences, tree depth, energy, $\widehat{R}$, and $\mathrm{ESS}$ all look good.
rstan::check_hmc_diagnostics(eval(fit)$fit)##
## Divergences:## 0 of 2000 iterations ended with a divergence.##
## Tree depth:## 0 of 2000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.if(max(rhat(eval(fit)), na.rm=T) >= 1.01) {
print("Warning: Rhat >=1.01")
} else {
print("All Rhat <1.01")
}## [1] "All Rhat <1.01"if(min(neff_ratio(eval(fit)), na.rm=T) <= 0.2) {
print("Warning: ESS <=0.2")
} else {
print("All ESS >0.2")
}## [1] "All ESS >0.2"2.3.9 Posterior Predictive Checks
We perform posterior predictive checks to assess how well the model fits the data.
pp_check(fit, type = "bars", nsamples = 250)## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.
2.3.10 Ploting
plot(fit)
## Extract and Plot Posterior Distributions
We extract posterior samples and plot the posterior distributions for the model coefficients.
posterior_samples <- posterior_samples(fit)## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.coefficients <- grep("b_", colnames(posterior_samples), value = TRUE)
plots <- lapply(coefficients, function(coef_name) {
ggplot(posterior_samples, aes_string(x = coef_name)) +
geom_density(fill = "blue", alpha = 0.5) +
labs(title = paste("Posterior Distribution for", gsub("b_", "", coef_name)),
x = "Effect Size", y = "Density")
})
#for (p in plots) {
# print(p)
#}2.3.11 Ploting the posterior probability
Next, let’s plot the posterior probability densities for our population-level estimates on the \(\logit\) scale.
parameter_names <- colnames(posterior_samples)
parameters <- grep("^b_", parameter_names, value = TRUE)
parameters <- parameters[!grepl("Intercept", parameters)]
print(parameters)## [1] "b_placementOptionM201_02" "b_placementOptionM201_03"
## [3] "b_placementOptionM201_04" "b_placementOptionM201_05"
## [5] "b_placementOptionM201_06" "b_placementOptionM201_07"
## [7] "b_placementOptionM201_08" "b_placementOptionM201_09"
## [9] "b_placementOptionM201_10" "b_ExperienceLevelLow"
## [11] "b_ExperienceLevelMedium" "b_ExpertiseEngineering"
## [13] "b_ExpertiseHumanFactors" "b_ExpertiseOther"
## [15] "b_AreaOfWorkNonMAutomotive" "b_AreaOfWorkOEM"
## [17] "b_AreaOfWorkResearchInstitute" "b_AreaOfWorkSupplier"
## [19] "b_Dev_StrucHybrid" "b_Dev_StrucWaterfall"parameters <- rev(parameters)
p <- mcmc_areas(fit, pars = parameters, prob_outer = 0.95, prob = 0.5)
#+ vline_0(size = 0.3)
y_labels <- c("Feature Definition", "Feature Realization", "System Design",
"Dedicated human factors team", "User-oriented teams", "Non-functional requirements team",
"System/feature evaluation team", "Safety team", "Person/team responsible for the overall system",
"ExperienceLevelLow", "ExperienceLevelMedium", "ExpertiseEngineering",
"ExpertiseHumanFactors", "ExpertiseOther", "AreaOfWorkNonMAutomotive",
"AreaOfWorkOEM", "AreaOfWorkResearchInstitute", "AreaOfWorkSupplier")
y_labels <- rev(y_labels)
p <- p + scale_y_discrete(labels = y_labels) ## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.print(p)
2.3.12 Conditional effects
Below we plot conditional effects for some of our estimates.
ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "placementOption", categorical = TRUE)
p_crit <- plot(p)[[1]] + # Removed plot = FALSE
scale_color_viridis(discrete = TRUE) +
scale_fill_viridis(discrete = TRUE) +
xlab("placementOption") +
ylab("") +
scale_x_continuous(breaks=seq(0,1))
ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Experience", categorical = TRUE)
p_crit <- plot(p)[[1]] + # Removed plot = FALSE
scale_color_viridis(discrete = TRUE) +
scale_fill_viridis(discrete = TRUE) +
xlab("Experience") +
ylab("") +
scale_x_continuous(breaks=seq(0,1))
ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Expertise", categorical = TRUE)
p_crit <- plot(p)[[1]] + # Removed plot = FALSE
scale_color_viridis(discrete = TRUE) +
scale_fill_viridis(discrete = TRUE) +
xlab("Criticality") +
ylab("") +
scale_x_continuous(breaks=seq(0,1))
3 Model Comparison
Given that our outcome variable follows an ordinal categorical structure, we have numerous modeling options available. Our outcome variable resembles that of an independent achievement model, akin to placing human factors expertise at different positions within an organization. Taking this characteristic into consideration, we anticipate that the following types of models may be suitable options for our analysis (Bürkner and Vuorre, 2019):
- Cumulative
- Adjacent-category (acat)
- Stoping ratio (sratio)
- Continuation ratio (cratio)
# Fit the model with `cumulative`
fit <- brm(
formula = model_formula,
data = modeldata,
family = cumulative,
prior = prior_spec,
chains = 4,
cores = 4,
seed = 12345
)## Warning: Rows containing NAs were excluded from the model.## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored## Start sampling## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 1 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -93.8941, but should be greater than the previous element, -93.8941 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 1## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 1 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -95.319, but should be greater than the previous element, -95.319 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 1## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 1 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -23.2631, but should be greater than the previous element, -23.2631 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 1## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 1 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 2 is -16679.4, but should be greater than the previous element, -16679.4 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 1## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 1 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -480.442, but should be greater than the previous element, -480.442 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 1## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 2 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -69.1526, but should be greater than the previous element, -69.1526 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 2## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 2 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -69.421, but should be greater than the previous element, -69.421 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 2## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 2 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 4 is -17.6883, but should be greater than the previous element, -17.6883 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 2## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 2 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is 6.00503e+258, but should be greater than the previous element, 6.00503e+258 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 2## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 2 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is 7.49579e+06, but should be greater than the previous element, 7.49579e+06 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 2## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 3 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -81.9327, but should be greater than the previous element, -81.9327 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 3## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 3 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -81.0321, but should be greater than the previous element, -81.0321 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 3## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 3 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 5 is -17.7521, but should be greater than the previous element, -17.7521 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 3## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 3 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 4 is 6.75391e+17, but should be greater than the previous element, 6.75391e+17 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 3## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 4 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 7 is -28.5643, but should be greater than the previous element, -28.5643 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 4## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 4 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 7 is -26.8515, but should be greater than the previous element, -26.8515 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 4## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 4 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is inf, but should be greater than the previous element, inf (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 4## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:## Chain 4 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 4 is 5.97504e+26, but should be greater than the previous element, 5.97504e+26 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.## Chain 4print(summary(fit))## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: ranking ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1 | ID)
## Data: modeldata (Number of observations: 256)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Group-Level Effects:
## ~ID (Number of levels: 30)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.17 0.14 0.01 0.53 1.00 2080 1952
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1] -3.00 0.62 -4.21 -1.81 1.00 1785
## Intercept[2] -2.12 0.60 -3.29 -0.95 1.00 1707
## Intercept[3] -1.48 0.59 -2.63 -0.35 1.00 1682
## Intercept[4] -0.95 0.59 -2.09 0.18 1.00 1706
## Intercept[5] -0.44 0.58 -1.59 0.67 1.00 1722
## Intercept[6] 0.07 0.58 -1.07 1.20 1.00 1777
## Intercept[7] 0.58 0.58 -0.55 1.71 1.00 1786
## Intercept[8] 1.18 0.58 0.05 2.33 1.00 1800
## Intercept[9] 2.16 0.60 1.02 3.34 1.00 1873
## placementOptionM201_02 -1.89 0.53 -2.90 -0.87 1.00 1898
## placementOptionM201_03 -0.61 0.52 -1.66 0.37 1.00 1867
## placementOptionM201_04 1.14 0.55 0.04 2.23 1.00 2178
## placementOptionM201_05 -0.98 0.53 -2.03 0.02 1.00 1899
## placementOptionM201_06 -1.47 0.51 -2.48 -0.46 1.00 1796
## placementOptionM201_07 1.05 0.54 -0.00 2.07 1.00 1900
## placementOptionM201_08 0.12 0.50 -0.87 1.11 1.00 1902
## placementOptionM201_09 -0.82 0.51 -1.81 0.17 1.00 1780
## placementOptionM201_10 -1.05 0.52 -2.05 -0.03 1.00 2005
## ExperienceLevelLow -0.32 0.51 -1.29 0.71 1.00 3582
## ExperienceLevelMedium -0.15 0.28 -0.69 0.41 1.00 3644
## ExpertiseEngineering -0.02 0.38 -0.75 0.73 1.00 2992
## ExpertiseHumanFactors 0.12 0.34 -0.56 0.78 1.00 3331
## ExpertiseOther -0.05 0.51 -1.07 0.97 1.00 4258
## AreaOfWorkNonMAutomotive -0.46 0.85 -2.20 1.14 1.00 2963
## AreaOfWorkOEM -0.11 0.41 -0.93 0.68 1.00 2175
## AreaOfWorkResearchInstitute 0.13 0.53 -0.91 1.19 1.00 2494
## AreaOfWorkSupplier 0.08 0.43 -0.78 0.93 1.00 2404
## Dev_StrucHybrid -0.10 0.36 -0.82 0.60 1.00 2895
## Dev_StrucWaterfall 0.08 0.36 -0.62 0.79 1.00 3133
## Tail_ESS
## Intercept[1] 2096
## Intercept[2] 2009
## Intercept[3] 2005
## Intercept[4] 2082
## Intercept[5] 1905
## Intercept[6] 2094
## Intercept[7] 2098
## Intercept[8] 2163
## Intercept[9] 2286
## placementOptionM201_02 2540
## placementOptionM201_03 2531
## placementOptionM201_04 2690
## placementOptionM201_05 2461
## placementOptionM201_06 2667
## placementOptionM201_07 2190
## placementOptionM201_08 2548
## placementOptionM201_09 2541
## placementOptionM201_10 2741
## ExperienceLevelLow 2891
## ExperienceLevelMedium 2947
## ExpertiseEngineering 2659
## ExpertiseHumanFactors 3052
## ExpertiseOther 2885
## AreaOfWorkNonMAutomotive 2656
## AreaOfWorkOEM 2163
## AreaOfWorkResearchInstitute 2605
## AreaOfWorkSupplier 2654
## Dev_StrucHybrid 2691
## Dev_StrucWaterfall 3017
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Fit the model with `cratio`
fit1 <- brm(
formula = model_formula,
data = modeldata,
family = cratio,
prior = prior_spec,
chains = 4,
cores = 4,
seed = 12345
)## Warning: Rows containing NAs were excluded from the model.## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored## Start samplingprint(summary(fit))## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: ranking ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1 | ID)
## Data: modeldata (Number of observations: 256)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Group-Level Effects:
## ~ID (Number of levels: 30)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.17 0.14 0.01 0.53 1.00 2080 1952
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1] -3.00 0.62 -4.21 -1.81 1.00 1785
## Intercept[2] -2.12 0.60 -3.29 -0.95 1.00 1707
## Intercept[3] -1.48 0.59 -2.63 -0.35 1.00 1682
## Intercept[4] -0.95 0.59 -2.09 0.18 1.00 1706
## Intercept[5] -0.44 0.58 -1.59 0.67 1.00 1722
## Intercept[6] 0.07 0.58 -1.07 1.20 1.00 1777
## Intercept[7] 0.58 0.58 -0.55 1.71 1.00 1786
## Intercept[8] 1.18 0.58 0.05 2.33 1.00 1800
## Intercept[9] 2.16 0.60 1.02 3.34 1.00 1873
## placementOptionM201_02 -1.89 0.53 -2.90 -0.87 1.00 1898
## placementOptionM201_03 -0.61 0.52 -1.66 0.37 1.00 1867
## placementOptionM201_04 1.14 0.55 0.04 2.23 1.00 2178
## placementOptionM201_05 -0.98 0.53 -2.03 0.02 1.00 1899
## placementOptionM201_06 -1.47 0.51 -2.48 -0.46 1.00 1796
## placementOptionM201_07 1.05 0.54 -0.00 2.07 1.00 1900
## placementOptionM201_08 0.12 0.50 -0.87 1.11 1.00 1902
## placementOptionM201_09 -0.82 0.51 -1.81 0.17 1.00 1780
## placementOptionM201_10 -1.05 0.52 -2.05 -0.03 1.00 2005
## ExperienceLevelLow -0.32 0.51 -1.29 0.71 1.00 3582
## ExperienceLevelMedium -0.15 0.28 -0.69 0.41 1.00 3644
## ExpertiseEngineering -0.02 0.38 -0.75 0.73 1.00 2992
## ExpertiseHumanFactors 0.12 0.34 -0.56 0.78 1.00 3331
## ExpertiseOther -0.05 0.51 -1.07 0.97 1.00 4258
## AreaOfWorkNonMAutomotive -0.46 0.85 -2.20 1.14 1.00 2963
## AreaOfWorkOEM -0.11 0.41 -0.93 0.68 1.00 2175
## AreaOfWorkResearchInstitute 0.13 0.53 -0.91 1.19 1.00 2494
## AreaOfWorkSupplier 0.08 0.43 -0.78 0.93 1.00 2404
## Dev_StrucHybrid -0.10 0.36 -0.82 0.60 1.00 2895
## Dev_StrucWaterfall 0.08 0.36 -0.62 0.79 1.00 3133
## Tail_ESS
## Intercept[1] 2096
## Intercept[2] 2009
## Intercept[3] 2005
## Intercept[4] 2082
## Intercept[5] 1905
## Intercept[6] 2094
## Intercept[7] 2098
## Intercept[8] 2163
## Intercept[9] 2286
## placementOptionM201_02 2540
## placementOptionM201_03 2531
## placementOptionM201_04 2690
## placementOptionM201_05 2461
## placementOptionM201_06 2667
## placementOptionM201_07 2190
## placementOptionM201_08 2548
## placementOptionM201_09 2541
## placementOptionM201_10 2741
## ExperienceLevelLow 2891
## ExperienceLevelMedium 2947
## ExpertiseEngineering 2659
## ExpertiseHumanFactors 3052
## ExpertiseOther 2885
## AreaOfWorkNonMAutomotive 2656
## AreaOfWorkOEM 2163
## AreaOfWorkResearchInstitute 2605
## AreaOfWorkSupplier 2654
## Dev_StrucHybrid 2691
## Dev_StrucWaterfall 3017
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Fit the model with `sratio`
fit2 <- brm(
formula = model_formula,
data = modeldata,
family = sratio,
prior = prior_spec,
chains = 4,
cores = 4,
seed = 12345
)## Warning: Rows containing NAs were excluded from the model.## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
##
## Start samplingprint(summary(fit))## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: ranking ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1 | ID)
## Data: modeldata (Number of observations: 256)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Group-Level Effects:
## ~ID (Number of levels: 30)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.17 0.14 0.01 0.53 1.00 2080 1952
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1] -3.00 0.62 -4.21 -1.81 1.00 1785
## Intercept[2] -2.12 0.60 -3.29 -0.95 1.00 1707
## Intercept[3] -1.48 0.59 -2.63 -0.35 1.00 1682
## Intercept[4] -0.95 0.59 -2.09 0.18 1.00 1706
## Intercept[5] -0.44 0.58 -1.59 0.67 1.00 1722
## Intercept[6] 0.07 0.58 -1.07 1.20 1.00 1777
## Intercept[7] 0.58 0.58 -0.55 1.71 1.00 1786
## Intercept[8] 1.18 0.58 0.05 2.33 1.00 1800
## Intercept[9] 2.16 0.60 1.02 3.34 1.00 1873
## placementOptionM201_02 -1.89 0.53 -2.90 -0.87 1.00 1898
## placementOptionM201_03 -0.61 0.52 -1.66 0.37 1.00 1867
## placementOptionM201_04 1.14 0.55 0.04 2.23 1.00 2178
## placementOptionM201_05 -0.98 0.53 -2.03 0.02 1.00 1899
## placementOptionM201_06 -1.47 0.51 -2.48 -0.46 1.00 1796
## placementOptionM201_07 1.05 0.54 -0.00 2.07 1.00 1900
## placementOptionM201_08 0.12 0.50 -0.87 1.11 1.00 1902
## placementOptionM201_09 -0.82 0.51 -1.81 0.17 1.00 1780
## placementOptionM201_10 -1.05 0.52 -2.05 -0.03 1.00 2005
## ExperienceLevelLow -0.32 0.51 -1.29 0.71 1.00 3582
## ExperienceLevelMedium -0.15 0.28 -0.69 0.41 1.00 3644
## ExpertiseEngineering -0.02 0.38 -0.75 0.73 1.00 2992
## ExpertiseHumanFactors 0.12 0.34 -0.56 0.78 1.00 3331
## ExpertiseOther -0.05 0.51 -1.07 0.97 1.00 4258
## AreaOfWorkNonMAutomotive -0.46 0.85 -2.20 1.14 1.00 2963
## AreaOfWorkOEM -0.11 0.41 -0.93 0.68 1.00 2175
## AreaOfWorkResearchInstitute 0.13 0.53 -0.91 1.19 1.00 2494
## AreaOfWorkSupplier 0.08 0.43 -0.78 0.93 1.00 2404
## Dev_StrucHybrid -0.10 0.36 -0.82 0.60 1.00 2895
## Dev_StrucWaterfall 0.08 0.36 -0.62 0.79 1.00 3133
## Tail_ESS
## Intercept[1] 2096
## Intercept[2] 2009
## Intercept[3] 2005
## Intercept[4] 2082
## Intercept[5] 1905
## Intercept[6] 2094
## Intercept[7] 2098
## Intercept[8] 2163
## Intercept[9] 2286
## placementOptionM201_02 2540
## placementOptionM201_03 2531
## placementOptionM201_04 2690
## placementOptionM201_05 2461
## placementOptionM201_06 2667
## placementOptionM201_07 2190
## placementOptionM201_08 2548
## placementOptionM201_09 2541
## placementOptionM201_10 2741
## ExperienceLevelLow 2891
## ExperienceLevelMedium 2947
## ExpertiseEngineering 2659
## ExpertiseHumanFactors 3052
## ExpertiseOther 2885
## AreaOfWorkNonMAutomotive 2656
## AreaOfWorkOEM 2163
## AreaOfWorkResearchInstitute 2605
## AreaOfWorkSupplier 2654
## Dev_StrucHybrid 2691
## Dev_StrucWaterfall 3017
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# Fit the model wuth `acat`
fit3 <- brm(
formula = model_formula,
data = modeldata,
family = acat,
prior = prior_spec,
chains = 4,
cores = 4,
seed = 12345
)## Warning: Rows containing NAs were excluded from the model.## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
##
## Start samplingprint(summary(fit))## Family: cumulative
## Links: mu = logit; disc = identity
## Formula: ranking ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1 | ID)
## Data: modeldata (Number of observations: 256)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Group-Level Effects:
## ~ID (Number of levels: 30)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.17 0.14 0.01 0.53 1.00 2080 1952
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1] -3.00 0.62 -4.21 -1.81 1.00 1785
## Intercept[2] -2.12 0.60 -3.29 -0.95 1.00 1707
## Intercept[3] -1.48 0.59 -2.63 -0.35 1.00 1682
## Intercept[4] -0.95 0.59 -2.09 0.18 1.00 1706
## Intercept[5] -0.44 0.58 -1.59 0.67 1.00 1722
## Intercept[6] 0.07 0.58 -1.07 1.20 1.00 1777
## Intercept[7] 0.58 0.58 -0.55 1.71 1.00 1786
## Intercept[8] 1.18 0.58 0.05 2.33 1.00 1800
## Intercept[9] 2.16 0.60 1.02 3.34 1.00 1873
## placementOptionM201_02 -1.89 0.53 -2.90 -0.87 1.00 1898
## placementOptionM201_03 -0.61 0.52 -1.66 0.37 1.00 1867
## placementOptionM201_04 1.14 0.55 0.04 2.23 1.00 2178
## placementOptionM201_05 -0.98 0.53 -2.03 0.02 1.00 1899
## placementOptionM201_06 -1.47 0.51 -2.48 -0.46 1.00 1796
## placementOptionM201_07 1.05 0.54 -0.00 2.07 1.00 1900
## placementOptionM201_08 0.12 0.50 -0.87 1.11 1.00 1902
## placementOptionM201_09 -0.82 0.51 -1.81 0.17 1.00 1780
## placementOptionM201_10 -1.05 0.52 -2.05 -0.03 1.00 2005
## ExperienceLevelLow -0.32 0.51 -1.29 0.71 1.00 3582
## ExperienceLevelMedium -0.15 0.28 -0.69 0.41 1.00 3644
## ExpertiseEngineering -0.02 0.38 -0.75 0.73 1.00 2992
## ExpertiseHumanFactors 0.12 0.34 -0.56 0.78 1.00 3331
## ExpertiseOther -0.05 0.51 -1.07 0.97 1.00 4258
## AreaOfWorkNonMAutomotive -0.46 0.85 -2.20 1.14 1.00 2963
## AreaOfWorkOEM -0.11 0.41 -0.93 0.68 1.00 2175
## AreaOfWorkResearchInstitute 0.13 0.53 -0.91 1.19 1.00 2494
## AreaOfWorkSupplier 0.08 0.43 -0.78 0.93 1.00 2404
## Dev_StrucHybrid -0.10 0.36 -0.82 0.60 1.00 2895
## Dev_StrucWaterfall 0.08 0.36 -0.62 0.79 1.00 3133
## Tail_ESS
## Intercept[1] 2096
## Intercept[2] 2009
## Intercept[3] 2005
## Intercept[4] 2082
## Intercept[5] 1905
## Intercept[6] 2094
## Intercept[7] 2098
## Intercept[8] 2163
## Intercept[9] 2286
## placementOptionM201_02 2540
## placementOptionM201_03 2531
## placementOptionM201_04 2690
## placementOptionM201_05 2461
## placementOptionM201_06 2667
## placementOptionM201_07 2190
## placementOptionM201_08 2548
## placementOptionM201_09 2541
## placementOptionM201_10 2741
## ExperienceLevelLow 2891
## ExperienceLevelMedium 2947
## ExpertiseEngineering 2659
## ExpertiseHumanFactors 3052
## ExpertiseOther 2885
## AreaOfWorkNonMAutomotive 2656
## AreaOfWorkOEM 2163
## AreaOfWorkResearchInstitute 2605
## AreaOfWorkSupplier 2654
## Dev_StrucHybrid 2691
## Dev_StrucWaterfall 3017
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc 1.00 0.00 1.00 1.00 NA NA NA
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Next, we’ll use Leave-One-Out Cross-Validation (LOO) to compare the models:
fit <- add_criterion(fit, criterion = "loo")
fit1 <- add_criterion(fit1, criterion = "loo")
fit2 <- add_criterion(fit2, criterion = "loo")
fit3 <- add_criterion(fit3, criterion = "loo")
loo_compare(fit, fit1, fit2, fit3)## elpd_diff se_diff
## fit1 0.0 0.0
## fit2 0.0 0.0
## fit -1.2 3.2
## fit3 -3.3 3.5The LOO analysis ranks “fit” as the top-performing model. Based on this ranking, it may be preferable to use “fit” for predictive modeling, as it demonstrates the best performance in terms of predictive accuracy according to the LOO comparison results. Note: We’ve demonstrated this process for “Effectiveness” here, but using the same script, you can perform it for other variables as well. We obtained similar results for all other variables.
4 Computational environment
print(sessionInfo(), locale=FALSE)## R version 4.3.3 (2024-02-29)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.4.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] ggthemes_5.1.0 viridis_0.6.5 viridisLite_0.4.2 patchwork_1.2.0
## [5] rethinking_2.31 cmdstanr_0.5.3 rstan_2.32.5 StanHeaders_2.32.6
## [9] mcmcplots_0.4.3 coda_0.19-4.1 mice_3.16.0 lubridate_1.9.3
## [13] forcats_1.0.0 stringr_1.5.1 purrr_1.0.2 readr_2.1.5
## [17] tidyr_1.3.1 tibble_3.2.1 tidyverse_2.0.0 brms_2.20.4
## [21] Rcpp_1.0.12 dplyr_1.1.4 ggplot2_3.5.0 gRain_1.4.1
## [25] gRbase_2.0.1 bnlearn_4.9.1 haven_2.5.4 bayesplot_1.11.1
##
## loaded via a namespace (and not attached):
## [1] tensorA_0.36.2.1 rstudioapi_0.15.0 jsonlite_1.8.8
## [4] shape_1.4.6.1 magrittr_2.0.3 jomo_2.7-6
## [7] farver_2.1.1 nloptr_2.0.3 rmarkdown_2.26
## [10] vctrs_0.6.5 minqa_1.2.6 denstrip_1.5.4
## [13] base64enc_0.1-3 htmltools_0.5.7 distributional_0.4.0
## [16] curl_5.2.1 broom_1.0.5 mitml_0.4-5
## [19] sass_0.4.8 bslib_0.6.1 htmlwidgets_1.6.4
## [22] plyr_1.8.9 zoo_1.8-12 cachem_1.0.8
## [25] sfsmisc_1.1-17 igraph_2.0.2 mime_0.12
## [28] lifecycle_1.0.4 iterators_1.0.14 pkgconfig_2.0.3
## [31] colourpicker_1.3.0 Matrix_1.6-5 R6_2.5.1
## [34] fastmap_1.1.1 shiny_1.8.0 digest_0.6.34
## [37] colorspace_2.1-0 ps_1.7.6 rprojroot_2.0.4
## [40] crosstalk_1.2.1 labeling_0.4.3 fansi_1.0.6
## [43] timechange_0.3.0 abind_1.4-5 compiler_4.3.3
## [46] withr_3.0.0 backports_1.4.1 inline_0.3.19
## [49] shinystan_2.6.0 highr_0.10 QuickJSR_1.1.3
## [52] pkgbuild_1.4.3 pan_1.9 MASS_7.3-60.0.1
## [55] gtools_3.9.5 loo_2.7.0 tools_4.3.3
## [58] httpuv_1.6.14 threejs_0.3.3 nnet_7.3-19
## [61] glue_1.7.0 nlme_3.1-164 promises_1.2.1
## [64] grid_4.3.3 checkmate_2.3.1 reshape2_1.4.4
## [67] generics_0.1.3 gtable_0.3.4 tzdb_0.4.0
## [70] data.table_1.15.2 hms_1.1.3 utf8_1.2.4
## [73] foreach_1.5.2 pillar_1.9.0 markdown_1.12
## [76] posterior_1.5.0 later_1.3.2 splines_4.3.3
## [79] lattice_0.22-5 survival_3.5-8 tidyselect_1.2.0
## [82] miniUI_0.1.1.1 knitr_1.45 gridExtra_2.3
## [85] V8_4.4.2 bookdown_0.39 stats4_4.3.3
## [88] xfun_0.42 bridgesampling_1.1-2 matrixStats_1.2.0
## [91] DT_0.32 stringi_1.8.3 yaml_2.3.8
## [94] boot_1.3-29 evaluate_0.23 codetools_0.2-19
## [97] cli_3.6.2 RcppParallel_5.1.7 rpart_4.1.23
## [100] shinythemes_1.2.0 xtable_1.8-4 processx_3.8.3
## [103] munsell_0.5.0 jquerylib_0.1.4 rstantools_2.4.0
## [106] ellipsis_0.3.2 dygraphs_1.1.1.6 Brobdingnag_1.2-9
## [109] lme4_1.1-35.1 glmnet_4.1-8 mvtnorm_1.2-4
## [112] ggridges_0.5.6 scales_1.3.0 xts_0.13.2
## [115] crayon_1.5.2 rlang_1.1.3 shinyjs_2.1.05 References
Bürkner, P.-C., Vuorre, M., 2019. Ordinal regression models in psychology: A tutorial. Advances in Methods and Practices in Psychological Science 2, 77–101. https://doi.org/10.1177/2515245918823199